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Abstract 

We describe a grid generation procedure designed to construct new classes 
of orthogonal coordinate systems for binary black hole spacetimes. The com- 
puted coordinates offer an alternative approach to current methods, in addi- 
tion to providing a framework for potentially more stable and accurate evolu- 
tions of colliding black holes. As a particular example, we apply our procedure 
to generate appropriate numerical grids to evolve Misner's axisymmetric ini- 
tial data set representing two equal mass black holes colliding head-on. These 
new results are compared with previously published calculations, and we find 
generally good agreement in both the waveform profiles and total radiated 
energies over the allowable range of separation parameters. Furthermore, be- 
cause no specialized treatment of the coordinate singularities is required, these 
new grids are more easily extendible to unequal mass and spinning black hole 
collisions. 
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I. INTRODUCTION 



After the first attempt in 1964 by Hahn and Lindquist [JTJ] to solve the problem of collid- 
ing black holes, a long-term and relatively successful program was initiated in the late 1960s 
by DeWitt et. al. (see DeWitt, Cadez, Smarr, and Eppley @-f)|], which we henceforth col- 
lectively refer to as DCSE). Their efforts were concentrated on simulating the axisymmetric 
head-on collision of two equal mass black holes using the time symmetric Misner data || for 
initial conditions. This data set possesses the "double Einstein-Rosen bridge" J7| topology 
in which two asymptotically flat sheets are joined by two throats representing non-rotating 
black holes. The throats are two-spheres that are invariant under an isometry operation 
identifying both sheets and, from a numerical standpoint, serve as a boundary across which a 
particular form of Neumann boundary condition must be imposed on the evolved quantities 
to preserve the isometry. 

DCSE developed the 2D Cadez |J coordinate system to numerically integrate and evolve 
the axisymmetric Einstein equations. The Cadez coordinates are curvilinear body-fitting 
coordinates that conform to the two black hole throats and become spherical at asymptotic 
infinity. The advantages of this coordinate system (and any other appropriate body-fitting 
system) include: (1) simplified boundary conditions at the black hole throats, (2) natural 
spherical-like grids which mimic the background Schwarzschild geometry that the solutions 
approach at late times, (3) simplification of waveform extraction calculations in the wave 
zone, and (4) exponential stretching of the grid in the radial direction which covers enough of 
the spacetime that asymptotic flatness can be applied to the conformal metric components at 
the outer grid edges. However, these conveniences are offset somewhat by the singular saddle 
point introduced by the coordinate transformation at the origin midway between the two 
black holes. DCSE attempted to deal with the saddle point by evolving the cylindrical (not 
Cadez) metric components everywhere on the Cadez grid and performing discretizations us- 
ing chain rule derivatives. Also, because black hole spacetimes can generate pathologically 
steep delta-function-like peaks in the metric components (assuming singularity avoiding 
time slices), numerical evolutions can quickly become unstable even if the inherent sym- 
metries are exploited in choosing the appropriate numerical grid and the evolved variables. 
Additionally, it is well known |||4| that an axis instability can be triggered in numerical 
evolutions of axisymmetric spacetimes. For all of these reasons, the computations of DCSE 
remained uncertain with error bars of order 100% ||. 

Recently, Anninos, Hobill, Seidel, Smarr and Suen l] (henceforth referred to as Pa- 
pers I, II and III respectively) improved upon the DCSE calculations by introducing a 
shift vector to diagonalize the Cadez metric and to evolve these components on the Cadez 
grid, thus suppressing the axis instability which is especially sensitive to the off-diagonal 
elements. The problem with the coordinate singularity at the origin was treated by con- 
structing a cylindrical coordinate "patch" to cover regions near the saddle point. The results 
of Anninos et al. are accurate to a few percent in the dominant wave signal for black holes 
which are initially placed at close to moderate separation distances. For larger initial sep- 
arations (greater than about 10M, where M is the single black hole mass), the evolutions 
become increasingly more inaccurate due in part to the fact that the saddle point remains 
within the causally connected regions of spacetime for longer periods of time as the black 
holes take longer to collide at the origin. If the black holes are initially separated by even 
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greater distances (more than about 20M), the evolutions become unstable and break down 
altogether. 

In short, the goal of long-term stable and accurate computations of highly separated 
colliding black holes has yet to be achieved. Here we propose two new classes of orthogonal 
body-fitting coordinate systems that are well-suited to the axisymmetric geometry of two 
colliding black holes, and that appear promising to improve upon current calculations. The 
new coordinates remove the singular saddle point from its obtrusive position at the origin 
and thus offer a potentially cleaner and more stable method of solution. Because we require 
spherical coordinate lines at the throats and in the far field, singular points are unavoidable. 
However, the saddle point singularities can be relocated along the z-axis to either the 
south or north pole of the top hole (and thus the north or south pole of the bottom hole, 
respectively), assuming the two holes are aligned vertically along the z-axis. We refer to 
the former system where the singularities face the opposite black hole as the class I system: 
class II refers to the latter case with the singularities facing away from the opposing throats. 

We are thus able to construct two new classes of coordinate systems, replacing the single 
singularity at the origin with a pair of saddle singularities on the throats. An advantage 
of transplanting the coordinate singularities to the throats prevents the black holes from 
"crashing" into a saddle point during the numerical evolutions. Instead, we implement a 
singularity avoiding lapse function that is zero on the throat, so that the evolution along the 
throat freezes in time. In addition, the (maximal) time slicing used in the evolutions will 
rapidly absorb coordinate lines into the event horizon and the saddle points will fall further 
inside the causally disconnected portion of the evolved spacetime where the lapse collapses 
exponentially to zero, providing an additional stabilizing element. 

The following sections describe our hybrid analytic/numerical procedure to generate dis- 
crete grids for axisymmetric binary black hole systems. Sec. [TT] introduces grid generating 
techniques in general and describes the particular analytic prescriptions that we developed 
for specifying one of the coordinates. Sec. JTJ describes the numerical calculation of the 
second coordinate orthogonal to the first, as well as the techniques used to compute the Ja- 
cobian matrix required for the coordinate transformations, and the tests which can be used 
to monitor the accuracy in which the coordinate systems and Jacobians are constructed. 
Results from actual numerical evolutions are presented in Sec. |IV| and compared with pre- 
viously published calculations. We summarize our results in Sec. [V]. 



II. GENERATING BODY FITTING COORDINATES: THE SPECIFIED 

COORDINATE 



A general and common method of constructing body-fitting coordinates is to let the 



curvilinear coordinates (77, £) satisfy the following elliptic partial differential equations [12 



d xx £ + d yy £ = P(r], f) , (la) 
d xx r] + d yy r] = Q(r], f) , (lb) 

where P(£,7/) and Q(C,,r]) are generating functions that can be used to adjust the behavior 
of the curvilinear coordinates in the physical domain. Because the boundaries are typically 
irregular in Cartesian coordinates (or cylindrical coordinates in our black hole work), and 
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the (77, £) coordinates are uniform in the transformed plane, it is desirable to carry out the 
computations in the transformed plane, switching the dependent and independent variables. 
Equations ([I]) can be inverted to yield 

adtfX — 2/3d V £X + 'yd^x + J 2 (Pd^x + Qd v x) = , (2a) 
ad^y - 2(3d ni y + jd vv y + J 2 (P% + Qd v y) = , (2b) 

where a = d^ + d^y 2 , (3 = d^xd v x+d v yd^y, 7 = d^x 2 + d^y 2 , and J = d^xd v y — d v xd^y. The 
problem of grid generation then reduces to solving a coupled set of nonlinear elliptic partial 
differential equations subject to the appropriate boundary conditions, which can themselves 
be intrinsically coupled in a complicated manner. Although the elliptic solution approach 
is a general and straight-forward one, other methods such as conformal mapping, algebraic 
transformations, and hyperbolic solutions of partial differential equations have also been 
developed (see |L2| for a review of grid generation techniques). 



In this paper, we present an altogether different and much simplified procedure that 
does not require one to solve coupled nonlinear elliptic equations. The idea rests on the 
notion that the grid generation process is greatly simplified if one of the coordinates can 
be specified analytically. Taking this coordinate as the one aligned asymptotically with 
either the radial or the angular direction, we show in this section how to construct a natural 
"specified" coordinate for three different classes of grids as characterized by the location of 
the singularity: class I (II) with two singular saddle points, one on each of the throats on the 
axis closest to (furthest from) the origin; and class III which are Cadez-like coordinates with 
a single singular point at the origin p = z = 0. We note that the function one chooses for 
the specified coordinate must solve the necessary boundary conditions, but is not restricted 
to satisfy any particular elliptic equation. 

Before we continue, a few notational comments are in order. Four axisymmetric coordi- 
nate systems are utilized in this paper: the combination (p, z) refers to standard cylindrical 
coordinates, (r, 6) denotes standard polar coordinates, (C £) are the radial- and angular- 
like body fitting coordinates, and (77, £) represent the logarithmic radial- and angular-like 
body fitting coordinates (note the two £ coordinates are identical). Finally we characterize 
the location of the two black holes by their radius a\ (02) and vertical distance from the 
origin along the z-axis Z\ (22), where Z\ (z 2 ) represents the absolute distance of the throat 
center from the origin, and the upper (lower) black hole is denoted by the subscript 1 (2). 



A. The Class I System 

The class I coordinate system has a coordinate (saddle-type) singularity where each of 
the throats meet the axis closest to the origin (see Fig. |l|). A singularity is generated at 
these points by requiring that a radial-like coordinate ( be zero on both of the throats and 
the entire section of the axis connecting the throats. Lines of constant ( will then transform 
from "peanut" -like surfaces near the throats to radial circles at infinity. The two throats 
and the axis between the throats in this case make up the line 77 = 0. The axis above (below) 
the top (bottom) hole is the angular coordinate value £ = (tt). In the equal mass case, the 
equator is the line £ = 9 = tt/2. The boundary conditions are shown in Fig. [| and Fig. [| in 
the (p,z) and (CO planes respectively. 
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FIG. 1. Class I coordinates for the n = 2.2 Misner data. 
An exponentially stretched radial coordinate satisfying the appropriate boundary condi- 
tions can be generated by first defining two radial distances from the two throats 



di = \Jp 2 + (z - zi) 2 - ai, 
d 2 



^ P 2 + {z + z 2 y 



0.2, 



(3) 
(4) 



and a third distance measure 
3 2 



p 2 + (z - z 1 + a l ) 2 + 



^p 2 + (z + z 2 - a 2 ) 2 - (zi +z 2 -0!- a 2 ) 



(5) 



defining an "elliptic" radial coordinate that is zero on the segment along the z-axis extending 
from z = z\ to z = — z 2 . (We note that the square root radicals implicitly refer to the absolute 
or positive root.) The actual curvilinear "radial" coordinates, ( and rj, are then constructed 
from these distance components by 



c 



(2/ci + l)did 2 d^ 
d\d 2 + Ki (did?, + d 2 dt) ' 



(6) 



and 



rj = sinh 1 I — 



K2 



(7) 




FIG. 2. The "radial" (Q) boundary conditions for the class I coordinate system. 
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FIG. 3. The class I schematic mapping of the rectangular (£, £) domain in relation to the 
throat, axis, and equator. The cross indicates the location of the coordinate singularity at the in- 
tersection of the throat and axis closest to the equator. 

Both are zero on the "spectacle" boundary of Fig. by construction. Equations (^) and (0) 
also have the proper behavior in the asymptotic limit where both p and z tend to infinity, 
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ie., ( — > r and rj — > log(r//«2)- The parameters K\ and /«2 are introduced so the coordinate 
systems can be "tuned" to improve the stability of evolutions. In particular, K\ controls 
the relative importance between the "elliptic" radius {d%) and the two distance-from-throat 
radii (d\ and d 2 ), and affects the shape of the grid close to the axis between the holes. By 
adjusting Ki, the grid can be pushed away from or drawn closer to the axis in this region, 
effectively increasing or decreasing the resolution between the holes in the p direction. k 2 
controls the asymptotic behavior of the coordinates in the far zone. Typical values for 
evolutions depend on the geometry of the system, and range from (k±, k 2 ) = (1, 2) for the 
/i = 1.2 case, to (0.7, 2) for the p = 2.7 case, where p is the Misner parameter defined in 
Sec. [TV^. 



B. The Class II System 

The class II coordinate system has a saddle point where each of the throats meet the 
axis furthest from the origin (see Fig. ^). A singularity is generated at these points by 
requiring that an angular-like coordinate £ be zero (tt) along the throat of the top (bottom) 
black hole and along the axis above (below) the top (bottom) throat. The coordinate £ also 
asymptotes to the polar angular coordinate 9. The north (south) throat in this case is a line 
of constant £ = (tc) between ( = and some finite ( value. The axis above (below) the 
north (south) throat is on the same constant £ line as the corresponding throats. The axis 
between the two throats is the line ( = 0. In the equal mass case, the equator is the line 
£ = 7r/2. The boundary conditions and mapping are shown in Fig. |5| and Fig. || in the (p, 
z) and (C, £) planes respectively 



2.5 r 




FIG. 4. Class II coordinates for the p = 2.2 Misner data. 
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FIG. 5. The "angular" (£) boundary conditions for the class II coordinates. 
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FIG. 6. The class II schematic mapping of the rectangular (CO domain in relation to the 
throat, axis, and equator. The cross indicates the location of the coordinate singularity at the in- 
tersection of the throat and axis furthest from the equator. 

The class II system can be generated by first defining two radial distances similar to the 
class I case 
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di a = p 2 + (z- Zl ) 2 - al, (8) 
d 2a = p 2 + {z + z 2 ) 2 - al (9) 



and two "hyperbolic" radial coordinates 



dib = -VP 2 + (z + z 2 + a 2 ) 



+ Jp 2 + (z - z 1 - ai) 2 + (zi + z 2 + ai + a 2 ), (10) 



d 2 b = +y P 2 + (z + z 2 + a 2 ) 2 



- yp 2 + (z - z 1 - ai) 2 + (zi + z 2 + a x + a 2 ), (11) 

which go to zero along the axis above and below the holes respectively. These coordinates 
can be combined to construct a function that is zero on the top hole's throat and along the 
axis above it (di), and a function that is zero on the bottom hole's throat and along the axis 
below it (d 2 ) 



* = + (12) 



, d 2a d 2b 

d 2 = 2 . , . — vj- (13) 

p 2 + (z + z 2 y 

The two functions d\ and d 2 are in turn combined to create a common function / that 
approaches +oo as d\ — > 0, and — oo as d 2 — > 

/ = (ai + a 2 + zi + z 2 ) Q- - — ) . (14) 

Finally, the function / is converted to an angular-like coordinate by 



^ArcCotp^fl + yTTT 5 ). (15) 

The angular coordinate in Eqn. (|l|) has the desired properties of being zero on the throat 
of the top hole and along the axis above it, tt on the bottom hole and along the axis below 
it, and approaches the normal spherical coordinate 9 in the asymptotic limit r — ► oo. 

In the class I system, the value of the second coordinate £ is specified on the outer grid 
edge to be the angular polar coordinate 9. There is slightly more freedom in class II for 
specifying the second coordinate rj at small r, since we only require that r = e v near the 
outer edge of the grid. We therefore define 

r = 2 sinh 77 + kii] (16) 

to be the radial coordinate. Because the throat is described by the first several 77 zones (see 
Fig. |j), k\ effectively controls the resolution along and near the throat boundary. 
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C. The Class III System: Cadez like Coordinates 



Cadez coordinates (see Fig. [j]) are related to cylindrical coordinates through the following 
semi-analytic complex transformation 



T) + it, = r [In (z + z + ip) + In (z - z + ip)] 



71=1 



+ 



(z + z + ip) n (z — z — ip) n 



(17) 



where z = ±zq are the locations of the throat centers. The constant radial and angular lines 
lie along the field and equipotential lines of two equally charged metallic cylinders located 
at the centers of the two throats. Hence the lines of constant r\ are spherical along the 
throats and in the asymptotic far field, and a singular saddle point is introduced midway 
between the two black holes at the origin p = z = 0. The coefficients C n are determined 
numerically by a least squares procedure to set the throats (defined by p\ h + (z t h ±-2o) 2 = a 2 , 
where a = a\ = ci2 is the throat radius) to lie on an rj = rjo = constant coordinate line. 




FIG. 7. Cadez coordinates for the Misner data with separation parameter fi = 2.2. 
It is relatively easy to generate a Cadez-like system using our grid generation prescrip- 
tion. A radial coordinate is required which vanishes on the throats (and only on the throats) 
and asymptotes to r as p and z tend to infinity. Clearly, this can be satisfied with 

did 2 



c 



where 



d\ + d 2 

d\ = 
d 2 = 



(l — k\ exp {^k 2 p 2 — k%z' 



p 2 + (z- z x ) 2 - Oi, 
p 2 + (z + z 2 ) 2 - a 2 , 



(18) 

(19) 
(20) 
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and the parameters k±, k2 and k% are introduced to control the behavior and resolution of 
the coordinates near the saddle point at the origin. In contrast with the Cadez coordinates, 
we note that the Cauchy-Riemann relations are not satisfied for this coordinate system. 



III. GENERATING THE SECOND COORDINATE AND OTHER NUMERICAL 

ISSUES 

We have developed two different methods for computing the second (unspecified) coor- 
dinate and to establish the positions of the (77, grid points in the cylindrical coordinate 
(p, z) space. In addition to evaluating the grid node positions, it is also necessary to com- 
pute the Jacobian matrix elements (ie., the coordinate derivatives) accurately. The Jacobian 
matrix and its inverse are needed to transform metric functions between the (p, z) and (77, 
£) spaces and to perform chain rule differentiations of the cylindrical metric functions on 
the curvilinear grids in the coordinate "patch" regions described in Sec. IVj. We turn first 
to the computation of the Jacobian, and then to the generation of the second coordinate. 



A. Computing the Jacobian 



In order to effectively use the curvilinear grids in numerical evolutions, it is important to 
generate the Jacobian matrix and its inverse as accurately as possible, since the gravitational 
waves emitted from the collisions are expected to be very weak in relation to the background 
geometry. Knowing the matrix components amounts to knowing, at each point on the 
grids, the coordinate derivatives dp/drj, dp/d£, dz/drj and dz/d^. A naive approach would 
be to simply evaluate these derivatives numerically at each point, after the grid has been 
calculated, using a second or higher order discretization stencil. Instead, we adopt a more 
accurate procedure and take advantage of the fact that the derivatives of the specified 
coordinate are known analytically. In this way, only two of the elements need to be computed 
numerically. 

The Jacobian matrix and its inverse can be written as 

d v p d v z 
d^p d^z 

where J = dpi] d z ^ — d z 7] <9 P £ is the Jacobian determinant. One can show from the orthogonal- 
ity of the (p, z) and (77, £ ) coordinates that the first derivatives obey the Cauchy-Riemann- 
like conditions 



1 
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-d z r] d p 7] 



(21) 



d z r] = adp£, (22a) 
d P v = -<rd& (22b) 



where a = y g m , and g^ is the metric in curvilinear coordinates. In the class I case, two 
components of the inverse Jacobian, d p rj and d z T), are known analytically. Thus, using the 
relationships (|2T| ) and ([22]), the following class I elements are determined exactly from the 
analytic expressions 
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U 71 

^' = iww' (23b) 

The remaining two elements, <9g;z and 9^p, are approximated by finite differencing the grid 
data using a centered fourth order stencil. In the class II case, the corresponding derivatives 
of the coordinate £ are known, and the expressions analogous to Eqs. ( |23| ) are derived 
simply by replacing rj with £. Additionally, in order to carry out the coordinate patch 
evolutions described in Sec. [IV A , it is also necessary to know the second derivatives of 



the coordinates. These are computed in the same manner as the Jacobian: the specified 
coordinate derivatives are evaluated analytically, and derivatives of the second coordinate 
are approximated numerically using fourth order stencils. 

We note that the orthogonality condition of the angular and radial coordinates 

or] ot, or] a£ 

can be imposed here to eliminate a third unknown derivative component from the Jacobian 
matrix, and reduce the number of numerical discretization operations to one. However, 
rather than taking this approach, we defer the use of the orthogonality relationship to test 
the accuracy in computing the coordinate systems and the Jacobian, and we have confirmed 
that Eqn. converges to zero with the order of the integration method. 



B. Line Integration in the (p, z) Plane 

One method to generate the second coordinate is to treat the lines of constant values as 
field lines of the known or specified coordinate. These field lines are computed by integrating 
along the normals to the iso-lines, e.g. along the gradient of the specified coordinate function. 
In other words, given a functional form for one of the coordinates, the second orthogonal 
coordinate, £ in the class I case, is determined by integrating the following equations 

(25a) 
(25b) 

along lines of constant £. Here A is an arbitrary integration parameter, and the normal- 
ization factor IV77I is introduced to keep the step sizes regular at points where 77 has large 
gradients. This procedure also follows for the class II coordinates, except that the integrated 
(differentiated) coordinate is r] (£). These lines (the second coordinate) are thus guaranteed 
to be orthogonal to the specified coordinate lines. 

Because the first coordinate is known exactly, we can form the spatial gradients ana- 
lytically, and the problem reduces conveniently to a straight-forward ODE integration. In 
practice, Eqs. ( ^5|) are solved using a 4th order Runge-Kutta integration with a small fixed 
step size. The only difficulty with this method comes with finding the (p, z) values at a 



dp 


d p r] 


dX 




dz 


d z 7] 


dX 
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particular (77, £) position, since the ODE integrations may overshoot the destination grid 
nodes. However, the step size is chosen small enough (typically about 1% of the width of 
a single grid zone) that we can interpolate linearly between the staggered points without 
sacrificing accuracy. Once these integrations are completed, the coordinate derivatives are 
evaluated using the techniques described in Sec. [Ill A| . 



C. Line Integration in the (77, £) Plane 



An alternate, and far more efficient, approach to constructing the numerical grids involves 
integrating the grid equations (0), or their class II equivalents 

diZ = oJ'+mr (26a) 
a ' p= oJimr (26b) 

in the (77, £) plane. That is, rather than supplying a starting (p, z) value, evaluating 77 (or 
£) and its gradients, and then tracing £ (or rf) along orthogonal lines in the (p, z) plane, we 
pre-suppose a regular (77, £) grid. For class I coordinates, we supply the (p, z) values along 
the outer boundary 77 = r] max which are consistent with a constant r surface in spherical 
coordinates, i.e., 

p = exp(77 max ) sin 9, (27a) 
z = exp(r] max ) cos 9, (27b) 



and £ = 9 which is evenly discretized along 77 = rj max . Eqs. (p3[) are then integrated inwards 
towards the throats. For the class II coordinates with equal mass black holes, we supply the 
(p, z) values along the equator and integrate Eqs. ( PE| ) towards the axis and throats. In this 
case, we use 

p = 2 sinliT/ + A477, (28a) 
z = 0, (28b) 

where 77 is evenly discretized from 77 = to r] max . Because the derivatives of the specified 
coordinates are known analytically, this integration scheme reduces to a set of ODEs that 
we solve using a second order Runge-Kutta scheme with a step size smaller or equal in size 
to the (77, £) grid spacing used in the evolutions. As a check on the accuracy of solutions, 
the resulting cylindrical coordinate values can be substituted into the analytic expression for 
fj(p,z) (class I) or £(p,z) (class II) to evaluate the accuracy of the numerical integrations. 
We find that they converge to the truncation order of the integration method. 



IV. APPLICATIONS 



In this final section we apply the class I grids to numerically solve the Einstein equations 
for the axisymmetric collision of two equal mass black holes, using the conformal Misner 
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solution for initial data. Although we have shown how to construct three different grid 
classes, we focus here only on the class I type: the class II evolutions have proven to be less 
stable than class I, and class III is similar to the Cadez case. To demonstrate the applicability 
and accuracy of the new class I grids to actual evolutions, we repeat the different parameter 
evolutions in Papers I — III and compare our results with the published calculations. 



A. Initial Data and Evolutions 



The Misner data set is an axisymmetric and time-symmetric, single parameter family of 
solutions with the conformally flat spatial 3-metric 



df 



dp 2 + dz 2 + p 2 



where p and z are the cylindrical coordinates, and 



\ sinh (n/i) 



+ 



1 



z + coth np) J p 2 + (z — coth np) x 



(29) 



(30) 



The conformal factor \l/ solves the Hamiltonian constraint with the proper isometry imposed 
between the upper and lower sheets and represents two equal mass, non-rotating black holes 
aligned along the axis of symmetry (the z-axis) , and centered at z = ± coth p with radius 
a = 1/ sinh p. The free parameter p defines the total or ADM mass of the spacetime and 
the proper distance along the spacelike geodesic connecting the two throats. Increasing p 
decreases the total mass of the system and sets the two holes further away from one another. 
For general axisymmetric transformed coordinates, the 3-metric (^9|) can be written as 



df 



c 



j(p 2 v + z 2 ) dr, 2 + j(pl + zf) di 2 + 



J sin 2 ^ 



(31) 



where # c = * /J l /\ J is a regularization variable that can be used to simplify the metric 
components and to provide a stabilizing element in the numerical evolutions, (rj, £) are 
the logarithmic radial-like and angular-like curvilinear coordinates, and the sin 2 £ term is 
explicitly factored out of the component partly for historical conventions and partly 
to help regularize the numerical evolutions: We evolve the conformal metric and extrinsic 
curvature components as described in Paper II. However, there J was defined to be the 
Jacobian determinant of the coordinate transformation {J = J) so that g m — — 1 
initially. Here we simply set J = 1 since there is no obvious advantage to renormalize 
one of the components when the Cauchy-Riemann conditions are not satisfied. The partial 
coordinate derivatives in ([H]) are computed as described in Sec. |III A| , completing the initial 
data. The data is then evolved according to the same procedures described in Paper II, 
ie., using the maximal slicing condition for the lapse function with antisymmetric boundary 
conditions at the throat surfaces, and an elliptic condition for the shift vector to preserve 
the metric in curvilinear coordinates to be diagonal throughout the evolution. 

To provide more stable evolutions, especially for the widely separated black hole cases, 
we utilize a coordinate "patch" as described in Paper II. Evolutions in this patched domain, 
which covers the saddle points and portions of the axis, are performed using the cylindrical 
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coordinate based metric components and chain rule derivatives to compute spatial gradients 
across the curvilinear grid nodes. The solutions are then transformed using the general tensor 
relations to reconstruct the curvilinear metric and extrinsic curvature components, and then 
linearly blending the results into the rest of the spacetime, which is evolved in a normal 
manner using the curvilinear metric components on the curvilinear grid. The coordinate 
patch used in Papers I — III, although only several zones deep in the angular direction, extends 
from the throat all the way out to the outer boundary in the radial direction. For comparison 
(see Fig. |8|), the domain in the class I coordinates which requires a patch is localized to 
just the first few zones in the radial direction. As a result, when the black holes merge 
to form a common event horizon, the patched coordinates and the associated numerical 
evolutions become irrelevant as the lapse collapses to zero over this region before the metric 
shear grows enough to disrupt the solutions. The evolutions with class I coordinates are 
therefore more robust and less sensitive to patch parameters then previous calculations. 




FIG. 8. Locations, shapes and sizes of the coordinate patches used in the Misner fi = 2.7 grids 
for the Cadez (horizontally hatched domain) and class I (diagonally hatched domain) cases. The 
shaded region labeled "radiation zone" represents the domain into which most of the gravitational 
radiation is emitted. That this zone is concentrated along the equator is attributed to the axisym- 
metric nature of the collision. Despite the representative shape of the radiation zone, the actual 
radiation extraction is performed on 2-spheres starting at a radius of 30M (= 15Madm) from the 
origin. 
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B. Results 
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FIG. 9. Gravitational waveforms emitted from the head-on collision of two black holes evolved 
with the Cadez and class I grids for the \x = 1.6, 2.2 and 3.0 Misner data. The waveforms 
are extracted at a radius of 30M from the origin, where M = Madm/^ is the approximate sin- 
gle black hole mass. The evolutions in each case are performed at a grid resolution of 200x 35 
(radialx angular) cells. 



We first consider three different cases with Misner parameters /x = 1.6, 2.2, and 3.0, and 
compare the class I gravitational waveforms with previous results using the Cadez grid. Each 
calculation is performed at the same grid resolution using 200x35 radialx angular zones. The 
Misner parameters are chosen so that the smaller value represents a data set in which the two 
black holes are already merged within a single event horizon, and the subsequent evolution 
is that of a single distorted black hole ringing down to the Schwarzschild solution as it emits 
gravitational waves. The data from the larger Misner parameter cases correspond to two 
distinct black holes separated by proper distances of L = 8.92M (/x = 2.2) and L = 15. 8M 
(/x = 3.0) between the two throats, where M is half the ADM mass of the spacetime (or 
approximately the mass of a single black hole). Comparison plots of the dominant £ = 2 
Zerilli waveforms, extracted at a distance of 30M from the origin, are shown in Fig. |9] for 
the three different Misner cases. We find maximum relative differences less than about 10% 
in amplitude and 4% in phase for the /x = 1.6 and 2.2 cases (with absolute deviations < 10 -2 
in amplitude), and up to about 200% and 4% differences in amplitude and phase for the 
more difficult /x = 3.0 case in which the black holes are initially highly separated. 

To emphasize the (in) stability of the solutions at late times when the signal crossing the 



detector becomes weaker, we plot in Fig. 10 the logarithm of the absolute value of the Zerilli 



function for the relatively uncertain /x = 3.0 case. Notice the class I grid solution maintains 
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a more regular oscillatory behavior and consistent damping rate throughout the wave signal 
and for longer periods of time than the Cadez case, which begins to break down at about 
70 M adm . 
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FIG. 10. The absolute value of the Zerilli function is shown for the u = 3.0 case using the 
class I and Cadez grids. The logarithmic scale highlights the oscillations in the waveform and the 
exponential damping of the wave amplitude. Although the waveforms compare favorably at smaller 
values of /i, there is a significant difference and improvement in the regularity of oscillations in the 
class I system compared with the less stable Cadez case for this large initial separation data. 




Next, we reproduce in Fig. 11 the equivalent of Fig. 14 in Paper III. The total £ = 2 



radiated energy (in units of the ADM mass, Madm = 2M) emitted during the black hole 
collisions is plotted as a function of the initial separation distance (in units of M) between 
the two throats. Also included in the figure are results from Paper III, the Davies et al. 
(DRPP) point particle calculation (E = 0.0104m 2 /M, plotted for m = M = M ADM /2), 
and its reduced mass correction (m — > mM/{m + M)). Results from the class I and Cadez 
grid evolutions match extremely well, better than 3%, for separation parameters /i < 2.2 
corresponding to initial physical separations of L < 8.92M, and improve significantly at 
smaller /i values. However, deviations between the class I and Cadez results become greater 
for black holes with further initial separations, reflecting the difficulty in evolving highly 
separated black holes for long periods of time. 

To assert a measure of uncertainty in our results, we perform several different calcula- 
tions of the /i = 2.7 and 3.0 data, varying the grid resolution, patch width, and coordinate 
parameters as defined in Sec. |II A| to manipulate the shape of the coordinate lines. For the 
more problematic Cadez evolutions, additional parameters include the patch length, dura- 
tion and diffusion, as well as varied treatments of the coordinate singularity (for example, 
shift vector specifications, discretization stencils, and regularization of certain metric and 
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curvature components). The symbols in Fig. |TT| represent the median results and the error 
bars indicate the variance with computational parameters. The variances are similar to the 
differences observed between the Cadez and class I results: roughly 30% and 100% in the 
[i = 2.7 and 3.0 cases respectively, with a trend for better agreement with increased resolu- 
tion (though the resolution studies are limited by the axis instability). For the [i < 2.2 cases 
the variances are too small to plot, but are consistent with the observed agreement in the 
evolutions. Considering the sensitivity of the results to grid and patch parameters, and the 
difficulties in evolving these systems using maximal slicing conditions and in axisymmetry, 
the overall results agree fairly well. Furthermore, the numerical calculations are in reason- 
able agreement with the reduced mass approximation and the semi-analytical calculation of 
Araiijo and Oliveira |14 in the large separation limit. 
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FIG. 11. Total gravitational wave energy (in units of the ADM mass) emitted from the head-on 
collision of two equal mass black holes as a function of the separation distance between the two 
throats (in units of M = Madm/2)- The class I and Cadez grid results are plotted together for 
comparison, along with the DRPP / [I3| / calculation of a particle falling into a more massive black 
hole, and its reduced mass correction. The error bars represent the uncertainties estimated by 
performing the evolutions for different computational parameters as described in the text. The un- 
certainties in both the class I and Cadez evolutions are comparable to differences between the two 
different grid results, ranging from < 3% for separations L < 10M, 30% for L ~ 13M , and 100% 
for L ~ 16M. 



Fig. [Ll] also demonstrates an added advantage of the new grid generation procedure 
to construct numerical grids in the very low Misner parameter cases. Due to convergence 
problems in the Newton-Raphson iterative inversion of Eqn. (|17D when the black hole 



throats are placed too close to the saddle point, previous calculations were limited to /i > 0.7 
| T5fl . The closest separation data shown in Fig. [ll] corresponds to /i = 0.5 (or L = 2.51M), 
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although grids for even smaller values of /i can be easily generated. 

An additional benefit from these new coordinates is their regularity at the origin, which 
makes calculations of the event horizon and null generators more accurate as the black holes 
merge. In Fig. |12| we show the evolution of the embedding of the event horizon found in 
the \i = 2.2 case evolved with the class I coordinates. The embedding of the horizon is 
smooth (except at the cusps on the z-axis), as were previous embeddings of the horizon 
in this spacetime [IB]. However, the null surfaces on the class I grid contain not only the 



horizon (in the domain 77 > 0), but also naturally contain the locus of generators waiting 
to find the horizon, as discussed in p7| . More precisely, to generate Fig. 11 of |T7] a 
coordinate transformation from Cadez to (p, z) coordinates was required. This coordinate 
transformation is not needed to evolve the locus in class I coordinates. This is more than a 
convenience, however. Since the null surface consisting of the horizon plus locus is naturally 
represented as a smooth continuous surface in Class I coordinates, the entire null surface 
can be unambiguously embedded in 3-space, which was not possible with the Cadez grid. 
This allows us to determine the separation between the horizons in embedded space using 
some relevant physical measure. While in Cadez coordinates the horizon separation was 
determined to keep the outer surface of the "pair of pants" figure smooth, here we determine 
the separation by embedding the entire null surface. This gives a natural separation between 
the horizons from the embedding of the locus and determines the geometry by the horizon. 
We use this separation to place the holes in the "wristwatch" (or "pair of pants from above" ) 
Fig. [12|. The class I coordinates also allow for more detailed examinations of the caustic line 
at the coalescence point, which will be discussed elsewhere. 
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FIG. 12. Embedding of the event horizon for the [i = 2.2 Misner data with the class I evolution. 
The geometry of the class I coordinates allows the separation between the holes in this figure to 
be physical, not artificial as in previous two black hole collision embedding diagrams derived from 
Cadez coordinates. Also, the regularity of the class I grid at the origin allows for more accurate 
examinations of the caustic line at the coalescence point. 
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V. CONCLUSIONS 



In many ways the axisymmetric evolution of two colliding black holes has been stalled 
due to the lack of a coordinate system without pathological instabilities or grid singularities. 
By developing new techniques to generate alternative grids, and by creating body-fitting 
grid geometries with singularities on the throats of the two black holes, we are able to 
achieve more stable long time evolutions of two black hole systems. We have demonstrated 
the applicability of these new grids in actual numerical evolutions of Misner's initial data 
set for the head-on collision of two equal mass black holes. The calculations presented here 
confirm the existing results to fairly good accuracy in the restricted stable parameter range, 
with maximum relative deviations less than 3% in the radiated energies for the fi < 2.2 
cases corresponding to initial proper separation distances of L < 8.92M (with significantly 
smaller deviations for the lower /i cases). The differences increase to about 30% for fi = 2.7 
(L = 12. 7M) and 100% for fi = 3.0 (L = 15. 8M), which are comparable to the uncertainties 
in the evolutions as defined in Sec. [IV B| for both the Cadez and class I systems. These 
differences reflect the difficulty in evolving black holes for long periods of time with maximal 
slicing conditions, and the sensitivity of the evolutions to treatments of the axis instability 
and coordinate singularities. 

Basically the Cadez and class I evolutions yield consistent results, even at the high /x val- 
ues, especially when considering the changes observed in waveforms and energies by varying 
the computational parameters. In addition to confirming previous results, and because no 
specialized treatment of the coordinate singularities is required, it seems promising that we 
can evolve physical systems with these new coordinates which were not as easily addressable 
with previous codes using either the Cadez or cylindrical coordinate systems. For example, 
we are currently extending this work to evolve spinning black hole and unequal mass black 
hole collisions. 

To conclude, we emphasize the following advantages of these new coordinates: (1) they 
achieve better zone coverage in the strong field interaction region near the origin - caus- 
tics, photon generators and embeddings of the event horizon are better resolved; (2) the 
coordinate patch extends over a much smaller domain in the class I system - just a few 
zones radially, and localized to the z-axis - so it is not an unstabilizing element at late 
times; (3) the resulting gravitational waveforms from evolutions are not as sensitive to the 
patch parameters, such as its width, length, duration and diffusion parameters; (4) due to 
its robustness and lack of a need for specialized treatments of the saddle points, the new 
code is more simplified and easily generalizable to include non-equal mass black holes and 
spinning black hole collisions; and (5) the new coordinates allow a larger range of initial data 
in the Misner parameter to be evolved, including evolutions of black holes that are further 
separated (though the accuracy is questionable for ji > 3.0), and more closely spaced (for 
smaller order perturbations) than in previous calculations. In addition, the new coordinates 
offer all the same advantages as the Cadez coordinates: They are logarithmic in the radial 
direction, and spherical on the throat and in the asymptotic wave zone, thus allowing for 
the same simplified treatment of boundary conditions and waveform extraction. 
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